import numpy as np

'''
此模块为各种变量之间的计算关系式
'''


# 水-水蒸气的相变潜热计算
def latent_heat(T):
    t = T - 273.15
    L = (2500 - 2.35 * t) * 1000
    return L


# 饱和水蒸气分压力计算（绝对温度）
def pressure_vap_sat(T):
    t = T - 273.15  # 转化为摄氏度
    Pv = 610.5 * np.exp(17.269 * t / (237.3 + t))  # 饱和水蒸气分压力计算公式
    return Pv


# 饱和水蒸气分压力关于温度的导数求解
def pressure_vap_sat_slope(T, flag=0, Pv_sat=None, L=None):
    if flag == 0:  # 差分方程求解导数
        h = 0.01
        delta = pressure_vap_sat(T + h) - pressure_vap_sat(T - h)
        Pv_sat_slope = delta / (2 * h)

    else:  # 克拉贝龙-克劳修斯方程求解导数
        Rv = 8.314472 / 0.018  # 水蒸气气体常数计算
        Pv_sat_slope = L * Pv_sat / (Rv * T ** 2)
    return Pv_sat_slope


# 材料的湿平衡曲线----------待完善
def U_curve(RH, a, b, c):
    U = a * RH / ((1 + b * RH) * (1 - c * RH))
    return U


# 等温吸放湿曲线斜率的计算
def U_slope(RH, a, b, c):
    h = 0.001  # 设置微元步长
    delta = U_curve((RH + h), a, b, c) - U_curve((RH - h), a, b, c)
    u_slope = delta / (2 * h)  # 利用差分计算导数
    return u_slope


# 系数D_v（水蒸气的渗透系数）求解------------李魁山拟合
def permeability_v(RH, d, e, f):
    D_v = d + e * (RH ** f)
    D_v /= 1000
    return D_v


# 系数D_l（液态水的渗透系数）的求解----------王莹莹算法
def permeability_l(T, RH, Pv_sat, D_v, den_l=1000):
    Rv = 8.314472 / 0.018  # 水蒸气气体常数计算
    D_l = (D_v * RH * Pv_sat) / (Rv * T * den_l)
    return D_l


# 水蒸气密度计算公式（温度，相对湿度）
def density_vapour(T, RH, Pv_sat):
    R = 8.314472  # 理想气体常数R约为8.31441 ± 0.00026 J/(mol·K)
    M = 0.018  # 水分子的摩尔质量为 18 g/mol，合0.018 kg/mol
    Rv = R / M  # 气态水分子的气体常数
    Pv = RH * Pv_sat  # 水蒸气分压力 = 相对湿度 * 饱和水蒸气分压力
    den_v = Pv / (Rv * T)  # 气体密度 = 气体分压力 / （气体常数 * 绝对温度）
    return den_v


# 系数ap_T（温度场时间项系数）求解
def ap_T(den_m, capacity_m):
    a_T = den_m * capacity_m  # 材料的总热容
    return a_T


# 系数bp_T（温度场扩散项系数）求解
def bp_T(T, RH, conductivity_eff, D_v):
    L = latent_heat(T)
    right = L * D_v * RH * pressure_vap_sat_slope(T)
    b_T = conductivity_eff + right  # bp_T 的计算公式
    return b_T


# 系数cp_T（温度耦合场扩散项系数）的求解
def cp_T(T, D_v, Pv_sat):
    L = latent_heat(T)
    c_T = L * D_v * Pv_sat  # cp_T 的计算公式
    return c_T


# 系数ap_W（湿度场时间项系数）求解
def ap_W(den_m, U_RH_slope):
    a_W = den_m * U_RH_slope
    return a_W


# 系数bp_W（湿度场扩散项系数）求解
def bp_W(T, RH, D_l, D_v, Pv_sat, den_l=1000):
    Rv = 8.314472 / 0.018  # 水蒸气气体常数计算
    c1 = D_v * Pv_sat
    c2 = D_l * den_l * Rv * T / RH
    b_W = c1 + c2
    return b_W


# 系数cp_W（湿度耦合场扩散项系数）求解
def cp_W(D_l, D_v, RH, Pv_sat_slope, den_l=1000):
    Rv = 8.314472 / 0.018  # 水蒸气气体常数计算
    c1 = D_l * den_l * Rv * np.log(RH)
    c2 = D_v * RH * Pv_sat_slope
    c_W = c1 + c2
    return c_W


# 温度场计算系数封装
def coefficients_T(T, RH,
                   den_m, conductivity_m, capacity_m,
                   d, e, f,
                   m_type):
    D_v = permeability_v(RH, d, e, f)  # 水蒸气扩散系数Dv计算
    Pv_sat = pressure_vap_sat(T)  # 饱和水蒸气分压力计算

    index = (m_type == 0)  # 找到空气层部分，空气层部分方程中Dv设置为0，以符合空气层传递公式
    if True in index:
        D_v[index] = 0

    a_T = ap_T(den_m, capacity_m)
    b_T = bp_T(T, RH, conductivity_m, D_v)
    c_T = cp_T(T, D_v, Pv_sat)
    return a_T, b_T, c_T


# 湿度场计算系数封装
def coefficients_W(T, RH, den_m, a, b, c, d, e, f, m_type):
    # 湿度场系数a计算
    U_RH_slope = U_slope(RH, a, b, c)  # 平衡曲线导数计算
    Pv_sat = pressure_vap_sat(T)  # 饱和水蒸气分压力Pv,s计算
    Rv = 8.314472 / 0.018  # 水蒸气气体常数Rv计算
    D_v = permeability_v(RH, d, e, f)  # 水蒸气扩散系数Dv计算
    D_l = permeability_l(T, RH, Pv_sat, D_v)  # 液态水扩散系数Dl计算
    Pv_sat_slope = pressure_vap_sat_slope(T)  # 水蒸气分压力导数d(Pv,s)计算
    # 湿度场系数a计算
    a_W = ap_W(den_m, U_RH_slope)  # 系数a计算

    index = (m_type == 0)  # 找到空气层部分
    if True in index:
        a_W[index] = Pv_sat[index] / (Rv * T[index])  # 空气层部分a计算
        D_l[index] = 0  # 空气层部分液态水传递系数设置为0

    # 湿度场系数b计算
    b_W = bp_W(T, RH, D_l, D_v, Pv_sat, den_l=1000)  # 系数b计算
    # 湿度场系数c计算
    c_W = cp_W(D_l, D_v, RH, Pv_sat_slope, den_l=1000)  # 系数c计算

    # 返回系数a,b,c
    return a_W, b_W, c_W


# 对流换热量
def q_convection(h, T_out, T_boundary):
    q = h * (T_out - T_boundary)
    return q


# 对流传质量
def j_convection(hm, T_out, T_boundary, W_out, W_boundary):
    Pv_sat_out = pressure_vap_sat(T_out)
    Pv_sat_boundary = pressure_vap_sat(T_boundary)
    den_v_out = density_vapour(T_out, W_out, Pv_sat_out)
    den_v_boundary = density_vapour(T_boundary, W_boundary, Pv_sat_boundary)
    j = hm * (den_v_out - den_v_boundary)
    return j


# 表面传质系数求解
def h_j(h):
    den_air = 1.20  # 空气密度
    cap_air = 1005  # 空气比热
    hm = h / (den_air * cap_air)  # 刘易斯关系求传质系数
    return hm


# 边界总换热量求解（无降雨情况）
def q_total(h, T_out, T_out_surface, W_out, W_out_surface, radiation=0, alpha=0.6):
    # 对流换热量求解
    qr = q_convection(h, T_out, T_out_surface)
    # 质量增量热量求解
    hm = h_j(h)
    ql = j_convection(hm, T_out, T_out_surface, W_out, W_out_surface) * latent_heat(T_out_surface)
    # 太阳辐射量求解
    qi = alpha * radiation
    # 总换热量求解：总换热量 = 对流换热量（qr） + 水蒸气内能增量(ql) + 太阳辐射量(qi)
    qt = qr + ql + qi
    return qt


# 表面换热系数h求解-----------待完善
def h_q(speed_wind=None):
    if speed_wind is None:
        h = 5
    else:
        h = 5 * speed_wind
    return h


# 相对湿度转化为水蒸气分压力
def rh2pressure(T, RH):
    P_v_s = pressure_vap_sat(T)
    P_v = RH * P_v_s
    return P_v


def main2():
    RH = 1
    a = np.array([0.055, 0.016, 0.082, 0.034, 0.088, 0.146])
    b = np.array([8.27, 4.55, 6.49, 6.69, 6.35, 0.814])
    c = np.array([0.63, 0.74, 0.645, 0.81, 0.53, 0.63])
    den_m = np.array([30, 42, 30, 1800, 2500, 1800])
    den_l = 1000

    mass_l = U_curve(RH, a, b, c)
    por_m = den_m * mass_l / den_l
    por_m = 1 / (1 / por_m + 1)
    print(por_m)


def main3():
    RH = np.linspace(0.2, 0.99, 100, dtype=float)
    d = 1.35e-10
    e = 5.46e-9
    f = 11.46
    a = permeability_v(RH, d, e, f)
    # T = np.linspace(273, 300, 100,dtype=float)
    # a, b, c = coefficients_W(T, RH, 440)

    import matplotlib.pyplot as plt
    plt.plot(RH, a)
    plt.show()


if __name__ == '__main__':
    print('*' * 20)
    den = 1000
    c = 4200
    a = ap_T(den, c)
    T = 89 + a
    print('*' * 20)
